! Project TWS Matrixexponentiële
!
!
! MODULE FFTW
!
! Deze module bevat functies om de ééndimensionale voorwaartse en inverse fouriertransformatie uit te voeren.
! De module gebruikt de library FFTW.
! 
! @author   Roel Matthysen
! @version  2.0


module fftw

    implicit none
    private

    integer, parameter :: doubleprecision = kind(0.d0)
    public :: fft,ifft
    
    contains

        ! Deze functie berekent de fouriertransformatie van een vector Ac
        ! 
        ! @param     Ac                    Een doubleprecision ééndimensionale matrix
        !
        ! @result    c                     De fouriergetransformeerde van Ac in een complex doubleprecision vector van
        !                                  dezelfde dimensie als Ac.
        !
        ! @since     2.0
        function fft(Ac) result(c)
            real(doubleprecision), dimension(:), intent(in) :: Ac
            complex(doubleprecision), dimension(size(Ac,1)) :: c
            integer*8 plan
            integer :: i,n
            n=size(Ac,1)
            call dfftw_plan_dft_r2c_1d(plan,n,Ac,c,'FFTW_ESTIMATE')
            call dfftw_execute_dft_r2c(plan, Ac, c)
            do i = 2,floor((size(Ac,1)+1)/2.0)
                c(n-i+2)=conjg(c(i))
            enddo
            call dfftw_destroy_plan(plan)
        
        end function
            
        ! Deze functie berekent de inverse fouriertransformatie van een vector Ac
        ! 
        ! @param     Ac                    Een complexe doubleprecision ééndimensionale matrix. Opdat het resultaat terug
        !                                  reëel is moet gelden dat Ac(i)=conjg(Ac(n-i))
        !
        ! @result    c                     De inverse fouriergetransformeerde van Ac in een reële doubleprecision vector
        !                                  van dezelfde dimensie als Ac.
        !
        ! @since     2.0
        function ifft(Ac) result(c)
            complex(doubleprecision), dimension(:), intent(in) :: Ac
            real(doubleprecision), dimension(size(Ac,1)) :: c
            integer*8 plan
            integer :: n
            n=size(Ac,1)
            call dfftw_plan_dft_c2r_1d(plan,n,Ac,c,'FFTW_ESTIMATE')
            call dfftw_execute_dft_c2r(plan, Ac, c)
            call dfftw_destroy_plan(plan)
            c=c/size(Ac,1)
            
        end function
    
end module
